home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / blas / dsymv.f < prev    next >
Text File  |  1997-06-25  |  8KB  |  266 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE DSYMV ( UPLO, N, ALPHA, A, LDA, X, INCX,
  5.      $                   BETA, Y, INCY )
  6. *     .. Scalar Arguments ..
  7.       DOUBLE PRECISION   ALPHA, BETA
  8.       INTEGER            INCX, INCY, LDA, N
  9.       CHARACTER*1        UPLO
  10. *     .. Array Arguments ..
  11.       DOUBLE PRECISION   A( LDA, * ), X( * ), Y( * )
  12. *     ..
  13. *
  14. *  Purpose
  15. *  =======
  16. *
  17. *  DSYMV  performs the matrix-vector  operation
  18. *
  19. *     y := alpha*A*x + beta*y,
  20. *
  21. *  where alpha and beta are scalars, x and y are n element vectors and
  22. *  A is an n by n symmetric matrix.
  23. *
  24. *  Parameters
  25. *  ==========
  26. *
  27. *  UPLO   - CHARACTER*1.
  28. *           On entry, UPLO specifies whether the upper or lower
  29. *           triangular part of the array A is to be referenced as
  30. *           follows:
  31. *
  32. *              UPLO = 'U' or 'u'   Only the upper triangular part of A
  33. *                                  is to be referenced.
  34. *
  35. *              UPLO = 'L' or 'l'   Only the lower triangular part of A
  36. *                                  is to be referenced.
  37. *
  38. *           Unchanged on exit.
  39. *
  40. *  N      - INTEGER.
  41. *           On entry, N specifies the order of the matrix A.
  42. *           N must be at least zero.
  43. *           Unchanged on exit.
  44. *
  45. *  ALPHA  - DOUBLE PRECISION.
  46. *           On entry, ALPHA specifies the scalar alpha.
  47. *           Unchanged on exit.
  48. *
  49. *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  50. *           Before entry with  UPLO = 'U' or 'u', the leading n by n
  51. *           upper triangular part of the array A must contain the upper
  52. *           triangular part of the symmetric matrix and the strictly
  53. *           lower triangular part of A is not referenced.
  54. *           Before entry with UPLO = 'L' or 'l', the leading n by n
  55. *           lower triangular part of the array A must contain the lower
  56. *           triangular part of the symmetric matrix and the strictly
  57. *           upper triangular part of A is not referenced.
  58. *           Unchanged on exit.
  59. *
  60. *  LDA    - INTEGER.
  61. *           On entry, LDA specifies the first dimension of A as declared
  62. *           in the calling (sub) program. LDA must be at least
  63. *           max( 1, n ).
  64. *           Unchanged on exit.
  65. *
  66. *  X      - DOUBLE PRECISION array of dimension at least
  67. *           ( 1 + ( n - 1 )*abs( INCX ) ).
  68. *           Before entry, the incremented array X must contain the n
  69. *           element vector x.
  70. *           Unchanged on exit.
  71. *
  72. *  INCX   - INTEGER.
  73. *           On entry, INCX specifies the increment for the elements of
  74. *           X. INCX must not be zero.
  75. *           Unchanged on exit.
  76. *
  77. *  BETA   - DOUBLE PRECISION.
  78. *           On entry, BETA specifies the scalar beta. When BETA is
  79. *           supplied as zero then Y need not be set on input.
  80. *           Unchanged on exit.
  81. *
  82. *  Y      - DOUBLE PRECISION array of dimension at least
  83. *           ( 1 + ( n - 1 )*abs( INCY ) ).
  84. *           Before entry, the incremented array Y must contain the n
  85. *           element vector y. On exit, Y is overwritten by the updated
  86. *           vector y.
  87. *
  88. *  INCY   - INTEGER.
  89. *           On entry, INCY specifies the increment for the elements of
  90. *           Y. INCY must not be zero.
  91. *           Unchanged on exit.
  92. *
  93. *
  94. *  Level 2 Blas routine.
  95. *
  96. *  -- Written on 22-October-1986.
  97. *     Jack Dongarra, Argonne National Lab.
  98. *     Jeremy Du Croz, Nag Central Office.
  99. *     Sven Hammarling, Nag Central Office.
  100. *     Richard Hanson, Sandia National Labs.
  101. *
  102. *
  103. *     .. Parameters ..
  104.       DOUBLE PRECISION   ONE         , ZERO
  105.       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  106. *     .. Local Scalars ..
  107.       DOUBLE PRECISION   TEMP1, TEMP2
  108.       INTEGER            I, INFO, IX, IY, J, JX, JY, KX, KY
  109. *     .. External Functions ..
  110.       LOGICAL            LSAME
  111.       EXTERNAL           LSAME
  112. *     .. External Subroutines ..
  113.       EXTERNAL           XERBLA
  114. *     .. Intrinsic Functions ..
  115.       INTRINSIC          MAX
  116. *     ..
  117. *     .. Executable Statements ..
  118. *
  119. *     Test the input parameters.
  120. *
  121.       INFO = 0
  122.       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
  123.      $         .NOT.LSAME( UPLO, 'L' )      )THEN
  124.          INFO = 1
  125.       ELSE IF( N.LT.0 )THEN
  126.          INFO = 2
  127.       ELSE IF( LDA.LT.MAX( 1, N ) )THEN
  128.          INFO = 5
  129.       ELSE IF( INCX.EQ.0 )THEN
  130.          INFO = 7
  131.       ELSE IF( INCY.EQ.0 )THEN
  132.          INFO = 10
  133.       END IF
  134.       IF( INFO.NE.0 )THEN
  135.          CALL XERBLA( 'DSYMV ', INFO )
  136.          RETURN
  137.       END IF
  138. *
  139. *     Quick return if possible.
  140. *
  141.       IF( ( N.EQ.0 ).OR.( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  142.      $   RETURN
  143. *
  144. *     Set up the start points in  X  and  Y.
  145. *
  146.       IF( INCX.GT.0 )THEN
  147.          KX = 1
  148.       ELSE
  149.          KX = 1 - ( N - 1 )*INCX
  150.       END IF
  151.       IF( INCY.GT.0 )THEN
  152.          KY = 1
  153.       ELSE
  154.          KY = 1 - ( N - 1 )*INCY
  155.       END IF
  156. *
  157. *     Start the operations. In this version the elements of A are
  158. *     accessed sequentially with one pass through the triangular part
  159. *     of A.
  160. *
  161. *     First form  y := beta*y.
  162. *
  163.       IF( BETA.NE.ONE )THEN
  164.          IF( INCY.EQ.1 )THEN
  165.             IF( BETA.EQ.ZERO )THEN
  166.                DO 10, I = 1, N
  167.                   Y( I ) = ZERO
  168.    10          CONTINUE
  169.             ELSE
  170.                DO 20, I = 1, N
  171.                   Y( I ) = BETA*Y( I )
  172.    20          CONTINUE
  173.             END IF
  174.          ELSE
  175.             IY = KY
  176.             IF( BETA.EQ.ZERO )THEN
  177.                DO 30, I = 1, N
  178.                   Y( IY ) = ZERO
  179.                   IY      = IY   + INCY
  180.    30          CONTINUE
  181.             ELSE
  182.                DO 40, I = 1, N
  183.                   Y( IY ) = BETA*Y( IY )
  184.                   IY      = IY           + INCY
  185.    40          CONTINUE
  186.             END IF
  187.          END IF
  188.       END IF
  189.       IF( ALPHA.EQ.ZERO )
  190.      $   RETURN
  191.       IF( LSAME( UPLO, 'U' ) )THEN
  192. *
  193. *        Form  y  when A is stored in upper triangle.
  194. *
  195.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  196.             DO 60, J = 1, N
  197.                TEMP1 = ALPHA*X( J )
  198.                TEMP2 = ZERO
  199.                DO 50, I = 1, J - 1
  200.                   Y( I ) = Y( I ) + TEMP1*A( I, J )
  201.                   TEMP2  = TEMP2  + A( I, J )*X( I )
  202.    50          CONTINUE
  203.                Y( J ) = Y( J ) + TEMP1*A( J, J ) + ALPHA*TEMP2
  204.    60       CONTINUE
  205.          ELSE
  206.             JX = KX
  207.             JY = KY
  208.             DO 80, J = 1, N
  209.                TEMP1 = ALPHA*X( JX )
  210.                TEMP2 = ZERO
  211.                IX    = KX
  212.                IY    = KY
  213.                DO 70, I = 1, J - 1
  214.                   Y( IY ) = Y( IY ) + TEMP1*A( I, J )
  215.                   TEMP2   = TEMP2   + A( I, J )*X( IX )
  216.                   IX      = IX      + INCX
  217.                   IY      = IY      + INCY
  218.    70          CONTINUE
  219.                Y( JY ) = Y( JY ) + TEMP1*A( J, J ) + ALPHA*TEMP2
  220.                JX      = JX      + INCX
  221.                JY      = JY      + INCY
  222.    80       CONTINUE
  223.          END IF
  224.       ELSE
  225. *
  226. *        Form  y  when A is stored in lower triangle.
  227. *
  228.          IF( ( INCX.EQ.1 ).AND.( INCY.EQ.1 ) )THEN
  229.             DO 100, J = 1, N
  230.                TEMP1  = ALPHA*X( J )
  231.                TEMP2  = ZERO
  232.                Y( J ) = Y( J )       + TEMP1*A( J, J )
  233.                DO 90, I = J + 1, N
  234.                   Y( I ) = Y( I ) + TEMP1*A( I, J )
  235.                   TEMP2  = TEMP2  + A( I, J )*X( I )
  236.    90          CONTINUE
  237.                Y( J ) = Y( J ) + ALPHA*TEMP2
  238.   100       CONTINUE
  239.          ELSE
  240.             JX = KX
  241.             JY = KY
  242.             DO 120, J = 1, N
  243.                TEMP1   = ALPHA*X( JX )
  244.                TEMP2   = ZERO
  245.                Y( JY ) = Y( JY )       + TEMP1*A( J, J )
  246.                IX      = JX
  247.                IY      = JY
  248.                DO 110, I = J + 1, N
  249.                   IX      = IX      + INCX
  250.                   IY      = IY      + INCY
  251.                   Y( IY ) = Y( IY ) + TEMP1*A( I, J )
  252.                   TEMP2   = TEMP2   + A( I, J )*X( IX )
  253.   110          CONTINUE
  254.                Y( JY ) = Y( JY ) + ALPHA*TEMP2
  255.                JX      = JX      + INCX
  256.                JY      = JY      + INCY
  257.   120       CONTINUE
  258.          END IF
  259.       END IF
  260. *
  261.       RETURN
  262. *
  263. *     End of DSYMV .
  264. *
  265.       END
  266.